This script explores whether two nuclear markers, HMGA2 and GATA6, are differentially expressed in KPC-T cells right next to acinar cells versus KPC-T cells distant from acinar cells, where no contact between the two cell types is possible.

In brief, murine acinar cells were isolated from wild-type C56Bl6/J mice and placed in cell culture in Waymouths’ medium. KPC-T cells were added to the acinar cells and the two cell types were allowed to coculture for a couple of days. Next, the cocultures were fixated and stained with GATA6 and HMGA2, with Hoechst as nuclear counterstain.

Images were acquired with Hoescht, HMGA2, tdTomato (endogenous for the KPC-T cells) and GATA6 signals. The image analysis in qupath consisted of the steps: 1. Annotation of all areas with KPC-T cells 2. Annotation of all acinar cells nearest to KPC-T cells (here called the “Interface”). In other words, this is the point where acinar cells and KPC-T cells meet. 3. Cell detection within the KPC-T area 4. Distance measure for all KPC-T cell objects to the Interface-annotation 5. Data export, with information for all individual KPC-T cells: the intensity of Hoechst, GATA6 and HMGA2, what well it came from, and the distance [µm] to the nearest point of acinar cells (nearest interface).

Of note, two filterings of the data are done in this code: - Two isolations were performed. However, the second isolation was later excluded from analysis since it was less successful. For example, the first isolation contained more acinar cells. - During imaging, we noticed that cells present in the middle of the culture wells had different expression profiles from the cells at the edges of the wells. The imaging settings were optimised to images first taken, which were in the middle of the wells. In general, the cells at the edges of the wells had a higher intensity (of everything), and often oversaturated. We conclude this likely stems from surface tension forming during incubation with antibodies in the wells rendered a higher coverage to the cells at the edges. Images from the middle of the wells are here called “Inner”, and images from the edges are called “outer”, and images from the “Outer” of wells are excluded.

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyr)
library(reshape2)
## 
## Attaching package: 'reshape2'
## The following object is masked from 'package:tidyr':
## 
##     smiths
library(ggthemes)
library(ggplot2)
library(rstatix)
## 
## Attaching package: 'rstatix'
## The following object is masked from 'package:stats':
## 
##     filter
library(ggpubr)
library(forcats)
setwd("/Users/sara.soderqvist/Library/CloudStorage/OneDrive-KarolinskaInstitutet/Writing_/Adaptive PDAC paper/Code to clean up and add in/read in") #Adjust working directory

coculture_og <- read.csv("read in_in vitro.csv")
coculture <- coculture_og
# In Qupath, where the image analysis was performed, some images had channel names with "..C2" (for HMGA2), "..C3" (for tdTomato) or "..C4" (for GATA6), and some did not have these suffixes. So same channel info ended up in different column for different cells. These are joined to one column here, the respective value in the proper channel column.

coculture_2 <- coculture %>%
   mutate(across(Nucleus..Hmga2.mean, ~ ifelse(is.na(.x), Nucleus..Hmga2..C2..mean, .x)))
coculture_3 <- coculture_2 %>%
   mutate(across(Nucleus..GATA6.mean, ~ ifelse(is.na(.x), Nucleus..GATA6..C4..mean, .x)))
coculture_4 <- coculture_3 %>%
   mutate(across(Cell..tdTomato.mean, ~ ifelse(is.na(.x), Cell..tdTomato..C3..mean, .x)))

na_test <- coculture_4[is.na(coculture_4$Cell..tdTomato.mean), ]
#na_test

#information in the column "Image":
#If it contains "control" = the cell comes from a well without acinar cells. 
#If it contains "intersec" or "interface" = the cell comes from a well with acinar cells, and from a field with an interfeace where KPC-T and acinar cells meet

coculture <- coculture_4[, -c(9:11)]
summary(is.na(coculture)) # There should only be True in Distance.. column. (from the "control" wells, which don't contain any acinar cells.) - ok
##    Image         Object.ID         Parent        Nucleus..Hoechst.mean
##  Mode :logical   Mode :logical   Mode :logical   Mode :logical        
##  FALSE:58510     FALSE:58510     FALSE:58510     FALSE:58510          
##                                                                       
##  Nucleus..Hmga2.mean Nucleus..GATA6.mean Cell..tdTomato.mean
##  Mode :logical       Mode :logical       Mode :logical      
##  FALSE:58510         FALSE:58510         FALSE:58510        
##                                                             
##  Distance.to.annotation.with.Interface.µm
##  Mode :logical                           
##  FALSE:51269                             
##  TRUE :7241

Normalising the nuclear HMGA2 and GATA6 expression to the to Hoechst signal for the individual cells:

coculture$HMGA2_fraction <-(coculture$Nucleus..Hmga2.mean/coculture$Nucleus..Hoechst.mean)
coculture$GATA6_fraction <-(coculture$Nucleus..GATA6.mean/coculture$Nucleus..Hoechst.mean)
coculture_cleanup <- separate(data = coculture, col = "Image", into="Image", sep=c(".nd2"), remove=TRUE)
## Warning: Expected 1 pieces. Additional pieces discarded in 58510 rows [1, 2, 3, 4, 5, 6,
## 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, ...].
coculture_cleanup$Setup <- (factor(nrow(coculture_cleanup)))

# Grouping of the KPCT cells based on if they are close to any acinar cells, or not in physical contact or near any acinar cell. The distance to the nearest acinar cell annotation is recorded in the "Distance.to.annotation.with.Interface.µm" column. The cutoff chosen for "near acinar cell annotation" is 20 µm.

coculture_cleanup$Setup <- "Above 20 µm"
coculture_cleanup <- coculture_cleanup %>% mutate(Setup = ifelse(coculture_cleanup$Distance.to.annotation.with.Interface.µm < 20, 'Within_20', .$Setup))

coculture_cleanup <- coculture_cleanup %>% mutate(Setup = ifelse(grepl('Control', coculture_cleanup$Parent), 'Ctrl_KPCT', .$Setup)) #Cells in this group are KPCT cultured in Waymouths medium only. They are included in some parts of the analysis, but not part of the primary research question included in the final plots for the manuscript.

#Some controls
na_test <- coculture_cleanup[is.na(coculture_cleanup$Setup), ]
is_not_na <- coculture_cleanup[!is.na(coculture_cleanup$Setup), ]
#na_test
#is_not_na
#unique(na_test$Image) 
coculture <- coculture_cleanup

Subgrouping based on the experimental setup, as mentioned above: Filter cells from images of the “inner” parts of the wells - this is the most accurate images. Images taken from the “outer” parts of the wells had a much stronger signal, moreover, these types of images were only captured from the second isolation.

coculture$Position <- (factor(nrow(coculture)))
coculture$Position <- "Inner"
coculture <- coculture %>% mutate(Position = ifelse(grepl('outer', coculture$Image), 'Outer', .$Position))

s_pos_coculture <- split(coculture, f = coculture$Position)
in_cocultures <- s_pos_coculture$Inner

Distribution of marker intensities stratified on the distance to acinar invasion front

#subset: only cells from "interface" images
interface_KPCT <- subset(coculture, Setup != "Ctrl_KPCT")
colnames(interface_KPCT) [8] <- "Distance_interface"
interface_KPCT <- interface_KPCT[, -c(3:7, 11, 12)]

melt_iKPCT <- melt(interface_KPCT, id = c("Image", "Object.ID", "Distance_interface"))

colnames(melt_iKPCT) [4] <- "Stain" 
colnames(melt_iKPCT) [5] <- "Cell_pos"

distrib_KPCTinterface <-ggplot(melt_iKPCT, aes(x=Distance_interface, y=Cell_pos, fill = Stain)) +
  geom_col(width = 0.7, position = "dodge", aes(alpha = 0.6)) +
  labs(title = "Distance to the acinar interface for each KPCT cell and its expression of HMGA2 and GATA6") +
  theme_clean() +
  facet_grid(Stain ~ .) +
  labs(x = "Distance to acinar interface [µm]", y = "Marker expression normalised to Hoechst") 

distrib_KPCTinterface
## Warning: `position_dodge()` requires non-overlapping x intervals.
## `position_dodge()` requires non-overlapping x intervals.

# see what is going on at the x = 0, where many cells seem to be (un-stack)
distrib_KPCTinterface_at30 <-ggplot(melt_iKPCT, aes(x=Distance_interface, y=Cell_pos, fill = Stain)) +
  geom_col(width = 0.7, position = "dodge", aes(alpha = 0.6)) +
  labs(title = "Marker expression first 30 µm") +
  theme_clean() +
  facet_grid(Stain ~ .) +
  labs(x = "Distance to acinar interface [µm]", y = "Marker expression normalised to Hoechst") +
  xlim(0, 30)

distrib_KPCTinterface_at30
## Warning: `position_dodge()` requires non-overlapping x intervals.
## `position_dodge()` requires non-overlapping x intervals.
## Warning: Removed 74086 rows containing missing values or values outside the scale range
## (`geom_col()`).

distrib_KPCTinterface_at100 <-ggplot(melt_iKPCT, aes(x=Distance_interface, y=Cell_pos, fill = Stain)) +
  geom_col(width = 0.7, position = "dodge", aes(alpha = 0.6)) +
  labs(title = "Marker expression first 100 µm") +
  theme_clean() +
  facet_grid(Stain ~ .) +
  labs(x = "Distance to acinar interface [µm]", y = "Marker expression normalised to Hoechst") +
  xlim(0, 100)

distrib_KPCTinterface_at100
## Warning: `position_dodge()` requires non-overlapping x intervals.
## Warning: `position_dodge()` requires non-overlapping x intervals.
## Warning: Removed 37452 rows containing missing values or values outside the scale range
## (`geom_col()`).

Display nr of cells in each group

s_in_cocultures <- split(in_cocultures, f = in_cocultures$Setup)
lapply(s_in_cocultures, nrow)
## $`Above 20 µm`
## [1] 39208
## 
## $Ctrl_KPCT
## [1] 4917
## 
## $Within_20
## [1] 9270
## Anova test

## Only Inner HMGA2
in_anova.test_hmga2 <- in_cocultures %>%
  anova_test(HMGA2_fraction ~ Setup) %>%
  adjust_pvalue(method = "BH") %>%
  add_significance("p.adj")
in_anova.test_hmga2
## ANOVA Table (type II tests)
## 
##   Effect DFn   DFd       F         p p<.05   ges     p.adj p.adj.signif
## 1  Setup   2 53392 595.092 2.47e-256     * 0.022 2.47e-256         ****
# tukeys as post-hoc

in_tukey_hmga2 <- in_cocultures %>% tukey_hsd(HMGA2_fraction~Setup) %>% add_xy_position()
in_tukey_hmga2
## # A tibble: 3 × 13
##   term  group1      group2    null.value estimate conf.low conf.high    p.adj
##   <chr> <chr>       <chr>          <dbl>    <dbl>    <dbl>     <dbl>    <dbl>
## 1 Setup Above 20 µm Ctrl_KPCT          0  -0.393   -0.420    -0.366  0       
## 2 Setup Above 20 µm Within_20          0  -0.0695  -0.0899   -0.0490 4.63e-14
## 3 Setup Ctrl_KPCT   Within_20          0   0.323    0.292     0.355  0       
## # ℹ 5 more variables: p.adj.signif <chr>, y.position <dbl>,
## #   groups <named list>, xmin <dbl>, xmax <dbl>
in_boxes_hmga2 <- ggboxplot(in_cocultures, x = "Setup", y = "HMGA2_fraction", 
          color = "Setup", fill = "Setup", alpha = 0.4, outlier.shape = NA) +
  xlab("Culture condition") + ylab("HMGA2 expression normalised to Hoechst signal") +
  stat_summary(fun.data = function(x) data.frame(y = 1.7, label = paste("Mean =",round(mean(x), digits = 4))), geom="text") +
 ylim(0, 2) +
  labs(title = "HMGA2 expression normalised to Hoechst signal - Inner images. Not displaying outliers")

in_boxes_hmga2_sign <- in_boxes_hmga2 + stat_pvalue_manual(label = "p.adj.signif",
                                                    in_tukey_hmga2, tip.length = 0.02)
in_boxes_hmga2_sign
## Warning: Removed 1998 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
## Warning: Removed 1998 rows containing non-finite outside the scale range
## (`stat_summary()`).
## Warning: Removed 3 rows containing non-finite outside the scale range
## (`stat_bracket()`).

## Only inner GATA6
in_anova.test_gata6 <- in_cocultures %>%
  anova_test(GATA6_fraction ~ Setup) %>%
  adjust_pvalue(method = "BH") %>%
  add_significance("p.adj")
in_anova.test_gata6
## ANOVA Table (type II tests)
## 
##   Effect DFn   DFd      F        p p<.05   ges    p.adj p.adj.signif
## 1  Setup   2 53392 86.955 1.98e-38     * 0.003 1.98e-38         ****
# tukeys as post-hoc
in_tukey_gata6 <- in_cocultures %>% tukey_hsd(GATA6_fraction~Setup) %>% add_xy_position()
in_tukey_gata6
## # A tibble: 3 × 13
##   term  group1      group2    null.value estimate conf.low conf.high       p.adj
##   <chr> <chr>       <chr>          <dbl>    <dbl>    <dbl>     <dbl>       <dbl>
## 1 Setup Above 20 µm Ctrl_KPCT          0   0.0272  0.0156     0.0388 0.000000125
## 2 Setup Above 20 µm Within_20          0   0.0479  0.0391     0.0568 0          
## 3 Setup Ctrl_KPCT   Within_20          0   0.0207  0.00719    0.0343 0.00097    
## # ℹ 5 more variables: p.adj.signif <chr>, y.position <dbl>,
## #   groups <named list>, xmin <dbl>, xmax <dbl>
in_boxes_gata6 <- ggboxplot(in_cocultures, x = "Setup", y = "GATA6_fraction", 
          color = "Setup", fill = "Setup", alpha = 0.4, outlier.shape = NA) +
  xlab("Culture condition") + ylab("GATA6 expression normalised to Hoechst signal") +
  stat_summary(fun.data = function(x) data.frame(y = 0.5, label = paste("Mean =",round(mean(x), digits = 4))), geom="text") +
   ylim(0, 0.7) +
  labs(title = "GATA6 expression normalised to Hoechst signal - Inner images. Not displaying outliers")

in_boxes_gata6_sign <- in_boxes_gata6 + stat_pvalue_manual(label = "p.adj.signif",
                                                    in_tukey_gata6, tip.length = 0.02)
in_boxes_gata6_sign
## Warning: Removed 2350 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
## Warning: Removed 2350 rows containing non-finite outside the scale range
## (`stat_summary()`).
## Warning: Removed 3 rows containing non-finite outside the scale range
## (`stat_bracket()`).

Additional test, using the absolute intensity values, i.e., not normalised to the Hoechst signal

## Only Inner HMGA2
in_abs_anova.test_hmga2 <- in_cocultures %>%
  anova_test(Nucleus..Hmga2.mean ~ Setup) %>%
  adjust_pvalue(method = "BH") %>%
  add_significance("p.adj")
in_abs_anova.test_hmga2
## ANOVA Table (type II tests)
## 
##   Effect DFn   DFd       F p p<.05   ges p.adj p.adj.signif
## 1  Setup   2 53392 844.025 0     * 0.031     0         ****
# tukeys as post-hoc

in_abs_tukey_hmga2 <- in_cocultures %>% tukey_hsd(Nucleus..Hmga2.mean~Setup) %>% add_xy_position()
in_abs_tukey_hmga2
## # A tibble: 3 × 13
##   term  group1      group2    null.value estimate conf.low conf.high      p.adj
##   <chr> <chr>       <chr>          <dbl>    <dbl>    <dbl>     <dbl>      <dbl>
## 1 Setup Above 20 µm Ctrl_KPCT          0   -163.   -173.      -153.  0         
## 2 Setup Above 20 µm Within_20          0     14.7     7.35      22.0 0.00000821
## 3 Setup Ctrl_KPCT   Within_20          0    178.    166.       189.  0         
## # ℹ 5 more variables: p.adj.signif <chr>, y.position <dbl>,
## #   groups <named list>, xmin <dbl>, xmax <dbl>
in_abs_boxes_hmga2 <- ggboxplot(in_cocultures, x = "Setup", y = "Nucleus..Hmga2.mean", 
          color = "Setup", fill = "Setup", alpha = 0.4, outlier.shape = NA) +
  xlab("Culture condition") + ylab("Absolute HMGA2 expression") +
  stat_summary(fun.data = function(x) data.frame(y = 800, label = paste("Mean =",round(mean(x), digits = 4))), geom="text") +
  ylim(0, 900) +
  labs(title = "Absolute HMGA2 expression - Inner images. Not displaying outliers")

in_abs_boxes_hmga2_sign <- in_abs_boxes_hmga2 + stat_pvalue_manual(label = "p.adj.signif",
                                                    in_abs_tukey_hmga2, tip.length = 0.02)
in_abs_boxes_hmga2_sign
## Warning: Removed 1342 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
## Warning: Removed 1342 rows containing non-finite outside the scale range
## (`stat_summary()`).
## Warning: Removed 3 rows containing non-finite outside the scale range
## (`stat_bracket()`).

## Only inner GATA6
in_abs_anova.test_gata6 <- in_cocultures %>%
  anova_test(Nucleus..GATA6.mean ~ Setup) %>%
  adjust_pvalue(method = "BH") %>%
  add_significance("p.adj")
in_abs_anova.test_gata6
## ANOVA Table (type II tests)
## 
##   Effect DFn   DFd       F         p p<.05   ges     p.adj p.adj.signif
## 1  Setup   2 53392 361.698 9.35e-157     * 0.013 9.35e-157         ****
# tukeys as post-hoc
in_abs_tukey_gata6 <- in_cocultures %>% tukey_hsd(Nucleus..GATA6.mean~Setup) %>% add_xy_position()
in_abs_tukey_gata6
## # A tibble: 3 × 13
##   term  group1      group2    null.value estimate conf.low conf.high    p.adj
##   <chr> <chr>       <chr>          <dbl>    <dbl>    <dbl>     <dbl>    <dbl>
## 1 Setup Above 20 µm Ctrl_KPCT          0     23.9     18.6      29.3 2.93e-14
## 2 Setup Above 20 µm Within_20          0     45.2     41.1      49.3 0       
## 3 Setup Ctrl_KPCT   Within_20          0     21.2     15.0      27.5 4.31e-14
## # ℹ 5 more variables: p.adj.signif <chr>, y.position <dbl>,
## #   groups <named list>, xmin <dbl>, xmax <dbl>
in_abs_boxes_gata6 <- ggboxplot(in_cocultures, x = "Setup", y = "Nucleus..GATA6.mean", 
          color = "Setup", fill = "Setup", alpha = 0.4, outlier.shape = NA) +
  xlab("Culture condition") + ylab("Absolute GATA6 expression") +
  stat_summary(fun.data = function(x) data.frame(y = 400, label = paste("Mean =",round(mean(x), digits = 4))), geom="text") +
  labs(title = "Absolute GATA6 expression - Inner images. Not displaying outliers") +
  ylim(0, 600)

in_abs_boxes_gata6_sign <- in_abs_boxes_gata6 + stat_pvalue_manual(label = "p.adj.signif",
                                                    in_abs_tukey_gata6, tip.length = 0.02)
in_abs_boxes_gata6_sign
## Warning: Removed 890 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
## Warning: Removed 890 rows containing non-finite outside the scale range
## (`stat_summary()`).
## Warning: Removed 3 rows containing non-finite outside the scale range
## (`stat_bracket()`).

Excluding the “Ctrl_KPCT” group, since it is not the most relevant for this question (they are non-coculture) Now the statistical test is a wilcoxon, since there are only two groups.

coculture_co <- subset(in_cocultures, Setup != "Ctrl_KPCT") 
unique(coculture_co$Setup)
## [1] "Above 20 µm" "Within_20"
coculture_co <- droplevels(coculture_co)
unique(coculture_co$Setup)
## [1] "Above 20 µm" "Within_20"
# Histograms of the groups' distributions 
gghistogram(coculture_co, x = "GATA6_fraction", color = "Setup", fill = "Setup", bins = 2000, palette = "Set1", title = "GATA6 histogram")

gghistogram(coculture_co, x = "HMGA2_fraction", color = "Setup", fill = "Setup", bins = 2000, palette = "Set1", title = "HMGA2 histogram")

## Only Inner HMGA2
in_wilcox_hmga2 <- coculture_co %>%
  wilcox_test(HMGA2_fraction ~ Setup) %>%
  adjust_pvalue(method = "BH") %>%
  add_significance("p.adj")

in_wilcox_hmga2 <- in_wilcox_hmga2 %>% add_xy_position()
in_wilcox_hmga2 [1, 10] <- 1.9
in_wilcox_hmga2
## # A tibble: 1 × 13
##   .y.         group1 group2    n1    n2 statistic        p    p.adj p.adj.signif
##   <chr>       <chr>  <chr>  <int> <int>     <dbl>    <dbl>    <dbl> <chr>       
## 1 HMGA2_frac… Above… Withi… 39208  9270 204987432 4.15e-82 4.15e-82 ****        
## # ℹ 4 more variables: y.position <dbl>, groups <named list>, xmin <dbl>,
## #   xmax <dbl>
in_co_boxes_hmga2 <- ggboxplot(coculture_co, x = "Setup", y = "HMGA2_fraction", 
          color = "Setup", fill = "Setup", alpha = 0.4, outlier.shape = NA) +
  xlab("Culture condition") + ylab("HMGA2 expression normalised to Hoechst signal") +
  stat_summary(fun.data = function(x) data.frame(y = 1.6, label = paste("Mean =",round(mean(x), digits = 4))), geom="text") +
 ylim(0, 2) +
  labs(title = "HMGA2 expression normalised to Hoechst signal - Inner images. Not displaying outliers") +
    scale_color_manual(values = c("#FC8D62", "#66C2A5")) +
  scale_fill_manual(values = c("#FC8D62", "#66C2A5"))

in_co_boxes_hmga2_sign <- in_co_boxes_hmga2 + stat_pvalue_manual(label = "p.adj.signif",
                                                    in_wilcox_hmga2, tip.length = 0.01)
in_co_boxes_hmga2_sign
## Warning: Removed 1987 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
## Warning: Removed 1987 rows containing non-finite outside the scale range
## (`stat_summary()`).

#As violin
in_co_violin_hmga2 <- ggviolin(coculture_co, x = "Setup", y = "HMGA2_fraction", 
          color = "Setup", fill = "Setup", alpha = 0.4, outlier.shape = NA) +
  xlab("Culture condition") + ylab("HMGA2 expression normalised to Hoechst signal") +
  stat_summary(fun.data = function(x) data.frame(y = 1.6, label = paste("Mean =",round(mean(x), digits = 4))), geom="text") +
 ylim(0, 2) +
  labs(title = "HMGA2 expression normalised to Hoechst signal - Inner images. Not displaying outliers") +
    scale_color_manual(values = c("#FC8D62", "#66C2A5")) +
  scale_fill_manual(values = c("#FC8D62", "#66C2A5"))

in_co_violin_hmga2_sign <- in_co_violin_hmga2 + stat_pvalue_manual(label = "p.adj.signif",
                                                    in_wilcox_hmga2, tip.length = 0.01)
in_co_violin_hmga2_sign
## Warning: Removed 1987 rows containing non-finite outside the scale range
## (`stat_ydensity()`).
## Removed 1987 rows containing non-finite outside the scale range
## (`stat_summary()`).
## Warning: Removed 121 rows containing missing values or values outside the scale range
## (`geom_violin()`).

## Only Inner GATA6
in_wilcox_gata6 <- coculture_co %>%
  wilcox_test(GATA6_fraction ~ Setup) %>%
  adjust_pvalue(method = "BH") %>%
  add_significance("p.adj")

in_wilcox_gata6 <- in_wilcox_gata6 %>% add_xy_position()
in_wilcox_gata6 [1, 10] <- 1.5
in_wilcox_gata6
## # A tibble: 1 × 13
##   .y.         group1 group2    n1    n2 statistic        p    p.adj p.adj.signif
##   <chr>       <chr>  <chr>  <int> <int>     <dbl>    <dbl>    <dbl> <chr>       
## 1 GATA6_frac… Above… Withi… 39208  9270 159932313 2.43e-72 2.43e-72 ****        
## # ℹ 4 more variables: y.position <dbl>, groups <named list>, xmin <dbl>,
## #   xmax <dbl>
in_co_boxes_gata6 <- ggboxplot(coculture_co, x = "Setup", y = "GATA6_fraction", 
          color = "Setup", fill = "Setup", alpha = 0.4, outlier.shape = NA) +
  xlab("Culture condition") + ylab("GATA6 expression normalised to Hoechst signal") +
  stat_summary(fun.data = function(x) data.frame(y = 0.9, label = paste("Mean =",round(mean(x), digits = 4))), geom="text") +
 ylim(0, 2) +
  labs(title = "GATA6 expression normalised to Hoechst signal - Inner images. Not displaying outliers") +
  scale_color_manual(values = c("#FC8D62", "#66C2A5")) +
  scale_fill_manual(values = c("#FC8D62", "#66C2A5"))

in_co_boxes_gata6_sign <- in_co_boxes_gata6 + stat_pvalue_manual(label = "p.adj.signif",
                                                    in_wilcox_gata6, tip.length = 0.01)
in_co_boxes_gata6_sign
## Warning: Removed 216 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
## Warning: Removed 216 rows containing non-finite outside the scale range
## (`stat_summary()`).

#As violin
in_co_violin_gata6 <- ggviolin(coculture_co, x = "Setup", y = "GATA6_fraction", 
          color = "Setup", fill = "Setup", alpha = 0.4, outlier.shape = NA) +
  xlab("Culture condition") + ylab("GATA6 expression normalised to Hoechst signal") +
  stat_summary(fun.data = function(x) data.frame(y = 0.9, label = paste("Mean =",round(mean(x), digits = 4))), geom="text") +
 ylim(0, 2) +
  labs(title = "GATA6 expression normalised to Hoechst signal - Inner images. Not displaying outliers") +
  scale_color_manual(values = c("#FC8D62", "#66C2A5")) +
  scale_fill_manual(values = c("#FC8D62", "#66C2A5"))

in_co_violin_gata6_sign <- in_co_violin_gata6 + stat_pvalue_manual(label = "p.adj.signif",
                                                    in_wilcox_gata6, tip.length = 0.01)
in_co_violin_gata6_sign
## Warning: Removed 216 rows containing non-finite outside the scale range
## (`stat_ydensity()`).
## Removed 216 rows containing non-finite outside the scale range
## (`stat_summary()`).
## Warning: Removed 42 rows containing missing values or values outside the scale range
## (`geom_violin()`).

s_in_co_cocultures <- split(coculture_co, f = coculture_co$Setup)
lapply(s_in_co_cocultures, nrow) 
## $`Above 20 µm`
## [1] 39208
## 
## $Within_20
## [1] 9270
# View how many wells that cells came from
# unique(s_in_co_cocultures[["Within_20"]]$Image)

Looking at distributions in individual wells and isolations:

coculture_co$Well_ID <- factor(nrow(coculture_co))
coculture_co <- coculture_co %>%
  mutate(Well_ID = ifelse(grepl('bckg_coculture plate well 1', .$Image), 'iso1_well1', .$Well_ID))
coculture_co <- coculture_co %>%
  mutate(Well_ID = ifelse(grepl('bckg_coculture plate well 2', .$Image), 'iso1_well2', .$Well_ID))
coculture_co <- coculture_co %>%
  mutate(Well_ID = ifelse(grepl('bckg_coculture plate well 3', .$Image), 'iso1_well3', .$Well_ID))
coculture_co <- coculture_co %>%
  mutate(Well_ID = ifelse(grepl('bckg_coculture plate well 4', .$Image), 'iso1_well4', .$Well_ID))
coculture_co <- coculture_co %>%
  mutate(Well_ID = ifelse(grepl('bckg_coculture plate well 5', .$Image), 'iso1_well5', .$Well_ID))
coculture_co <- coculture_co %>%
  mutate(Well_ID = ifelse(grepl('bckg_coculture plate 2_well 1', .$Image), 'iso2_well1', .$Well_ID))
coculture_co <- coculture_co %>%
  mutate(Well_ID = ifelse(grepl('bckg_coculture plate 2_well 2', .$Image), 'iso2_well2', .$Well_ID))
coculture_co <- coculture_co %>%
  mutate(Well_ID = ifelse(grepl('bckg_coculture plate 2_well 3', .$Image), 'iso2_well3', .$Well_ID))
coculture_co <- coculture_co %>%
  mutate(Well_ID = ifelse(grepl('bckg_coculture plate 2_well 4', .$Image), 'iso2_well4', .$Well_ID))


gghistogram(coculture_co, x = "GATA6_fraction", color = "Well_ID", fill = "Well_ID", facet.by = "Setup", bins = 2000, palette = "Set1", title = "GATA6 histogram colored on well ID")

gghistogram(coculture_co, x = "HMGA2_fraction", color = "Well_ID", fill = "Well_ID", facet.by = "Setup", bins = 2000, palette = "Set1", title = "HMGA2 histogram colored on well ID")

ggdotplot(coculture_co, y = "GATA6_fraction", x = "Setup", color = "Well_ID", fill = "Well_ID", binwidth = 0.1, palette = "Set1", title = "GATA6 dotplot colored on well ID", label.rectangle = TRUE)

ggdotplot(coculture_co, y = "HMGA2_fraction", x = "Setup", color = "Well_ID", fill = "Well_ID", binwidth = 0.1, palette = "Set1", title = "HMGA2 dotplot colored on well ID", label.rectangle = TRUE)

## Only Inner GATA6, grouped by the well
in_wilcox_gata6_g <- coculture_co %>%
  group_by(Well_ID) %>%
  wilcox_test(GATA6_fraction ~ Setup) %>%
  adjust_pvalue(method = "BH") %>%
  add_significance("p.adj")

in_wilcox_gata6_g <- in_wilcox_gata6_g %>% add_xy_position()
in_wilcox_gata6_g[c(1:9), 11] <- 1.5
in_wilcox_gata6_g
## # A tibble: 9 × 14
##   Well_ID    .y.           group1 group2    n1    n2 statistic        p    p.adj
##   <chr>      <chr>         <chr>  <chr>  <int> <int>     <dbl>    <dbl>    <dbl>
## 1 iso1_well1 GATA6_fracti… Above… Withi…  4521   446    879827 8.89e- 6 2.00e- 5
## 2 iso1_well2 GATA6_fracti… Above… Withi…  6002  1215   3619262 6.84e- 1 6.84e- 1
## 3 iso1_well3 GATA6_fracti… Above… Withi…  9566  2033   9644035 5.61e- 1 6.31e- 1
## 4 iso1_well4 GATA6_fracti… Above… Withi…  7525  1181   4517786 3.55e- 1 4.56e- 1
## 5 iso1_well5 GATA6_fracti… Above… Withi…  6884  1844   5661927 1.01e-12 9.09e-12
## 6 iso2_well1 GATA6_fracti… Above… Withi…  1028   682    363003 2.13e- 1 3.20e- 1
## 7 iso2_well2 GATA6_fracti… Above… Withi…  1735   926    891198 3.23e- 6 1.45e- 5
## 8 iso2_well3 GATA6_fracti… Above… Withi…  1108   443    281799 4.98e- 6 1.49e- 5
## 9 iso2_well4 GATA6_fracti… Above… Withi…   839   500    233850 4.3 e- 4 7.74e- 4
## # ℹ 5 more variables: p.adj.signif <chr>, y.position <dbl>,
## #   groups <named list>, xmin <dbl>, xmax <dbl>
in_co_boxes_gata6_g <- ggboxplot(coculture_co, x = "Setup", y = "GATA6_fraction", 
          color = "Setup", fill = "Setup", alpha = 0.4, outlier.shape = NA, facet.by = "Well_ID") +
  xlab("Culture condition") + ylab("GATA6 expression normalised to Hoechst signal") +
  stat_summary(fun.data = function(x) data.frame(y = 0.9, label = paste("Median =",round(median(x), digits = 4))), geom="text") +
 ylim(0, 2) +
  labs(title = "GATA6 expression normalised to Hoechst signal - Inner images. Not displaying outliers", subtitle = "unadjusted p values per well") +
  scale_color_manual(values = c("#FC8D62", "#66C2A5")) +
  scale_fill_manual(values = c("#FC8D62", "#66C2A5"))

in_co_boxes_gata6_g_sign <- in_co_boxes_gata6_g + stat_pvalue_manual(label = "p",
                                                    in_wilcox_gata6_g, tip.length = 0.01)
in_co_boxes_gata6_g_sign
## Warning: Removed 216 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
## Warning: Removed 216 rows containing non-finite outside the scale range
## (`stat_summary()`).

## Only Inner HMGA2, grouped by the well
in_wilcox_HMGA2_g <- coculture_co %>%
  group_by(Well_ID) %>%
  wilcox_test(HMGA2_fraction ~ Setup) %>%
  adjust_pvalue(method = "BH") %>%
  add_significance("p.adj")

in_wilcox_HMGA2_g <- in_wilcox_HMGA2_g %>% add_xy_position()
in_wilcox_HMGA2_g[c(1:9), 11] <- 3
in_wilcox_HMGA2_g
## # A tibble: 9 × 14
##   Well_ID    .y.            group1 group2    n1    n2 statistic        p   p.adj
##   <chr>      <chr>          <chr>  <chr>  <int> <int>     <dbl>    <dbl>   <dbl>
## 1 iso1_well1 HMGA2_fraction Above… Withi…  4521   446   1171467 1.59e- 8 7.15e-8
## 2 iso1_well2 HMGA2_fraction Above… Withi…  6002  1215   3923418 2.85e- 5 6.41e-5
## 3 iso1_well3 HMGA2_fraction Above… Withi…  9566  2033  10583956 3.54e-10 3.19e-9
## 4 iso1_well4 HMGA2_fraction Above… Withi…  7525  1181   4419650 7.66e- 1 7.66e-1
## 5 iso1_well5 HMGA2_fraction Above… Withi…  6884  1844   6393020 6.32e- 1 7.11e-1
## 6 iso2_well1 HMGA2_fraction Above… Withi…  1028   682    363743 1.87e- 1 2.40e-1
## 7 iso2_well2 HMGA2_fraction Above… Withi…  1735   926    898950 4.06e- 7 1.22e-6
## 8 iso2_well3 HMGA2_fraction Above… Withi…  1108   443    234261 1.61e- 1 2.40e-1
## 9 iso2_well4 HMGA2_fraction Above… Withi…   839   500    197191 6.65e- 2 1.20e-1
## # ℹ 5 more variables: p.adj.signif <chr>, y.position <dbl>,
## #   groups <named list>, xmin <dbl>, xmax <dbl>
in_co_boxes_HMGA2_g <- ggboxplot(coculture_co, x = "Setup", y = "HMGA2_fraction", 
          color = "Setup", fill = "Setup", alpha = 0.4, outlier.shape = NA, facet.by = "Well_ID") +
  xlab("Culture condition") + ylab("HMGA2 expression normalised to Hoechst signal") +
  stat_summary(fun.data = function(x) data.frame(y = 2.5, label = paste("Median =",round(median(x), digits = 4))), geom="text") +
 ylim(0, 5) +
  labs(title = "HMGA2 expression normalised to Hoechst signal - Inner images. Not displaying outliers", subtitle = "unadjusted p values per well") +
  scale_color_manual(values = c("#FC8D62", "#66C2A5")) +
  scale_fill_manual(values = c("#FC8D62", "#66C2A5"))

in_co_boxes_HMGA2_g_sign <- in_co_boxes_HMGA2_g + stat_pvalue_manual(label = "p",
                                                    in_wilcox_HMGA2_g, tip.length = 0.01)
in_co_boxes_HMGA2_g_sign
## Warning: Removed 229 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
## Warning: Removed 229 rows containing non-finite outside the scale range
## (`stat_summary()`).

Plots with medians instead of means displayed:

in_co_vioin_gata6_g <- ggviolin(coculture_co, x = "Setup", y = "GATA6_fraction", 
          color = "Setup", fill = "Setup", alpha = 0.4, outlier.shape = NA, facet.by = "Well_ID") +
  xlab("Culture condition") + ylab("GATA6 expression normalised to Hoechst signal") +
  stat_summary(fun.data = function(x) data.frame(y = 0.9, label = paste("Median =",round(median(x), digits = 4))), geom="text") +
 ylim(0, 2) +
  labs(title = "GATA6 expression normalised to Hoechst signal - Inner images. Not displaying outliers", subtitle = "unadjusted p values per well") +
  scale_color_manual(values = c("#FC8D62", "#66C2A5")) +
  scale_fill_manual(values = c("#FC8D62", "#66C2A5"))

in_co_vioin_gata6_g_sign <- in_co_vioin_gata6_g + stat_pvalue_manual(label = "p",
                                                    in_wilcox_gata6_g, tip.length = 0.01)
in_co_vioin_gata6_g_sign
## Warning: Removed 216 rows containing non-finite outside the scale range
## (`stat_ydensity()`).
## Warning: Removed 216 rows containing non-finite outside the scale range
## (`stat_summary()`).
## Warning: Removed 711 rows containing missing values or values outside the scale range
## (`geom_violin()`).

in_co_violin_HMGA2_g <- ggviolin(coculture_co, x = "Setup", y = "HMGA2_fraction", 
          color = "Setup", fill = "Setup", alpha = 0.4, outlier.shape = NA, facet.by = "Well_ID") +
  xlab("Culture condition") + ylab("HMGA2 expression normalised to Hoechst signal") +
  stat_summary(fun.data = function(x) data.frame(y = 2.5, label = paste("Median =",round(median(x), digits = 4))), geom="text") +
 ylim(0, 5) +
  labs(title = "HMGA2 expression normalised to Hoechst signal - Inner images. Not displaying outliers", subtitle = "unadjusted p values per well") +
  scale_color_manual(values = c("#FC8D62", "#66C2A5")) +
  scale_fill_manual(values = c("#FC8D62", "#66C2A5"))

in_co_violin_HMGA2_g_sign <- in_co_violin_HMGA2_g + stat_pvalue_manual(label = "p",
                                                    in_wilcox_HMGA2_g, tip.length = 0.01)
in_co_violin_HMGA2_g_sign
## Warning: Removed 229 rows containing non-finite outside the scale range
## (`stat_ydensity()`).
## Warning: Removed 229 rows containing non-finite outside the scale range
## (`stat_summary()`).
## Warning: Removed 461 rows containing missing values or values outside the scale range
## (`geom_violin()`).

in_co_vioin_gata6_g <- ggviolin(coculture_co, x = "Setup", y = "GATA6_fraction", 
          color = "Setup", fill = "Setup", alpha = 0.4, outlier.shape = NA, facet.by = "Well_ID") +
  xlab("Culture condition") + ylab("GATA6 expression normalised to Hoechst signal") +
  stat_summary(fun.data = function(x) data.frame(y = 0.9, label = paste("Mean =",round(mean(x), digits = 4))), geom="text") +
 ylim(0, 2) +
  labs(title = "GATA6 expression normalised to Hoechst signal - Inner images. Not displaying outliers", subtitle = "unadjusted p values per well") +
  scale_color_manual(values = c("#FC8D62", "#66C2A5")) +
  scale_fill_manual(values = c("#FC8D62", "#66C2A5"))

in_co_vioin_gata6_g_sign <- in_co_vioin_gata6_g + stat_pvalue_manual(label = "p",
                                                    in_wilcox_gata6_g, tip.length = 0.01)
in_co_vioin_gata6_g_sign
## Warning: Removed 216 rows containing non-finite outside the scale range
## (`stat_ydensity()`).
## Warning: Removed 216 rows containing non-finite outside the scale range
## (`stat_summary()`).
## Warning: Removed 711 rows containing missing values or values outside the scale range
## (`geom_violin()`).

in_co_violin_HMGA2_g <- ggviolin(coculture_co, x = "Setup", y = "HMGA2_fraction", 
          color = "Setup", fill = "Setup", alpha = 0.4, outlier.shape = NA, facet.by = "Well_ID") +
  xlab("Culture condition") + ylab("HMGA2 expression normalised to Hoechst signal") +
  stat_summary(fun.data = function(x) data.frame(y = 2.5, label = paste("Mean =",round(mean(x), digits = 4))), geom="text") +
 ylim(0, 5) +
  labs(title = "HMGA2 expression normalised to Hoechst signal - Inner images. Not displaying outliers", subtitle = "unadjusted p values per well") +
  scale_color_manual(values = c("#FC8D62", "#66C2A5")) +
  scale_fill_manual(values = c("#FC8D62", "#66C2A5"))

in_co_violin_HMGA2_g_sign <- in_co_violin_HMGA2_g + stat_pvalue_manual(label = "p",
                                                    in_wilcox_HMGA2_g, tip.length = 0.01)
in_co_violin_HMGA2_g_sign
## Warning: Removed 229 rows containing non-finite outside the scale range
## (`stat_ydensity()`).
## Warning: Removed 229 rows containing non-finite outside the scale range
## (`stat_summary()`).
## Warning: Removed 461 rows containing missing values or values outside the scale range
## (`geom_violin()`).

Plots grouped on the isolation rounds

coculture_co$Isolation <- factor(nrow(coculture_co))
coculture_co <- coculture_co %>%
  mutate(Isolation = ifelse(grepl('iso1', .$Well_ID), 'iso_1', .$Isolation))
coculture_co <- coculture_co %>%
  mutate(Isolation = ifelse(grepl('iso2', .$Well_ID), 'iso_2', .$Isolation))
unique(coculture_co$Isolation)
## [1] "iso_1" "iso_2"
gghistogram(coculture_co, x = "GATA6_fraction", color = "Isolation", fill = "Isolation", facet.by = "Setup", bins = 2000, palette = "Set1", title = "GATA6 histogram colored on Isolation")

gghistogram(coculture_co, x = "HMGA2_fraction", color = "Isolation", fill = "Isolation", facet.by = "Setup", bins = 2000, palette = "Set1", title = "HMGA2 histogram colored on Isolation")

ggdotplot(coculture_co, y = "GATA6_fraction", x = "Setup", color = "Isolation", fill = "Isolation", binwidth = 0.1, palette = "Set1", title = "GATA6 dotplot colored on Isolation", label.rectangle = TRUE)

ggdotplot(coculture_co, y = "HMGA2_fraction", x = "Setup", color = "Isolation", fill = "Isolation", binwidth = 0.1, palette = "Set1", title = "HMGA2 dotplot colored on Isolation", label.rectangle = TRUE)

## Only Inner GATA6, grouped by the Isolation
in_wilcox_gata6_iso <- coculture_co %>%
  group_by(Isolation) %>%
  wilcox_test(GATA6_fraction ~ Setup) %>%
  adjust_pvalue(method = "BH") %>%
  add_significance("p.adj")

in_wilcox_gata6_iso <- in_wilcox_gata6_iso %>% add_xy_position()
in_wilcox_gata6_iso[c(1:2), 11] <- 2
in_wilcox_gata6_iso
## # A tibble: 2 × 14
##   Isolation .y.            group1   group2    n1    n2 statistic       p   p.adj
##   <chr>     <chr>          <chr>    <chr>  <int> <int>     <dbl>   <dbl>   <dbl>
## 1 iso_1     GATA6_fraction Above 2… Withi… 34498  6719 111939188 9.23e-6 9.23e-6
## 2 iso_2     GATA6_fraction Above 2… Withi…  4710  2551   6424876 9.91e-7 1.98e-6
## # ℹ 5 more variables: p.adj.signif <chr>, y.position <dbl>,
## #   groups <named list>, xmin <dbl>, xmax <dbl>
in_co_boxes_gata6_iso <- ggboxplot(coculture_co, x = "Setup", y = "GATA6_fraction", 
          color = "Setup", fill = "Setup", alpha = 0.4, outlier.shape = NA, facet.by = "Isolation") +
  xlab("Culture condition") + ylab("GATA6 expression normalised to Hoechst signal") +
  stat_summary(fun.data = function(x) data.frame(y = 0.9, label = paste("Median =",round(median(x), digits = 4))), geom="text") +
 ylim(0, 2.5) +
  labs(title = "GATA6 expression normalised to Hoechst signal - Inner images. Not displaying outliers", subtitle = "unadjusted p values per well") +
  scale_color_manual(values = c("#FC8D62", "#66C2A5")) +
  scale_fill_manual(values = c("#FC8D62", "#66C2A5"))

in_co_boxes_gata6_iso_sign <- in_co_boxes_gata6_iso + stat_pvalue_manual(label = "p",
                                                    in_wilcox_gata6_iso, tip.length = 0.01)
in_co_boxes_gata6_iso_sign
## Warning: Removed 118 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
## Warning: Removed 118 rows containing non-finite outside the scale range
## (`stat_summary()`).

## Only Inner HMGA2, grouped by the Isolation
in_wilcox_HMGA2_iso <- coculture_co %>%
  group_by(Isolation) %>%
  wilcox_test(HMGA2_fraction ~ Setup) %>%
  adjust_pvalue(method = "BH") %>%
  add_significance("p.adj")

in_wilcox_HMGA2_iso <- in_wilcox_HMGA2_iso %>% add_xy_position()
in_wilcox_HMGA2_iso[c(1:2), 11] <- 4
in_wilcox_HMGA2_iso
## # A tibble: 2 × 14
##   Isolation .y.            group1   group2    n1    n2 statistic       p   p.adj
##   <chr>     <chr>          <chr>    <chr>  <int> <int>     <dbl>   <dbl>   <dbl>
## 1 iso_1     HMGA2_fraction Above 2… Withi… 34498  6719 121209302 2.61e-9 5.22e-9
## 2 iso_2     HMGA2_fraction Above 2… Withi…  4710  2551   6084151 3.69e-1 3.69e-1
## # ℹ 5 more variables: p.adj.signif <chr>, y.position <dbl>,
## #   groups <named list>, xmin <dbl>, xmax <dbl>
in_co_boxes_HMGA2_iso <- ggboxplot(coculture_co, x = "Setup", y = "HMGA2_fraction", 
          color = "Setup", fill = "Setup", alpha = 0.4, outlier.shape = NA, facet.by = "Isolation") +
  xlab("Culture condition") + ylab("HMGA2 expression normalised to Hoechst signal") +
  stat_summary(fun.data = function(x) data.frame(y = 2.5, label = paste("Median =",round(median(x), digits = 4))), geom="text") +
 ylim(0, 5) +
  labs(title = "HMGA2 expression normalised to Hoechst signal - Inner images. Not displaying outliers", subtitle = "unadjusted p values per well") +
  scale_color_manual(values = c("#FC8D62", "#66C2A5")) +
  scale_fill_manual(values = c("#FC8D62", "#66C2A5"))

in_co_boxes_HMGA2_iso_sign <- in_co_boxes_HMGA2_iso + stat_pvalue_manual(label = "p",
                                                    in_wilcox_HMGA2_iso, tip.length = 0.01)
in_co_boxes_HMGA2_iso_sign
## Warning: Removed 229 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
## Warning: Removed 229 rows containing non-finite outside the scale range
## (`stat_summary()`).

Using a mean value from each well to create boxplots:

  1. get means from every well (per: Marker (HMGA2 and GATA6), Setup (= within or above 20 µm), well_ID) to a well-based data frame
  2. run test and make boxplot with these as the new groups
HMGA2og <- coculture_co[, c(2, 9, 11, 13, 14)]
GATA6og <- coculture_co[, c(2, 10, 11, 13, 14)]
## HMGA2
#Calculate means of each well

s_HMGA2og <- split(HMGA2og, f = HMGA2og$Well_ID)

means_HMGA2og <- lapply(X = s_HMGA2og, FUN = function (x) {
  split_df <- split(x, f = x$Setup)
 # head(split_df)
  lapply(X = split_df, FUN = function (a) {
    h_mean <- median(a$HMGA2_fraction)
    h_mean
  })
}) #This function now prints out all means for each well_ID and Setup.
#means_HMGA2og
#Example for one of the wells:
means_HMGA2og$iso1_well1
## $`Above 20 µm`
## [1] 0.4387586
## 
## $Within_20
## [1] 0.3414553
means_HMGA2bound <- bind_rows(means_HMGA2og, .id = "Well_ID")

means_HMGA2bm <- melt(means_HMGA2bound, id.vars = "Well_ID", value.name = "HMGA2_mean", variable.name = "Setup")

means_HMGA2bm <- means_HMGA2bm[order(means_HMGA2bm[,1] ),] # Order to enable a paired test

## GATA6
s_GATA6og <- split(GATA6og, f = GATA6og$Well_ID)

means_GATA6og <- lapply(X = s_GATA6og, FUN = function (x) {
  split_df <- split(x, f = x$Setup)
 # head(split_df)
  lapply(X = split_df, FUN = function (a) {
    h_mean <- median(a$GATA6_fraction)
    h_mean
  })
}) #This function now prints out all means for each well_ID and Setup.
#means_GATA6og
#Example for one of the wells:
means_GATA6og$iso1_well1
## $`Above 20 µm`
## [1] 0.1149244
## 
## $Within_20
## [1] 0.128859
means_GATA6bound <- bind_rows(means_GATA6og, .id = "Well_ID")

means_GATA6bm <- melt(means_GATA6bound, id.vars = "Well_ID", value.name = "GATA6_mean", variable.name = "Setup")

means_GATA6bm <- means_GATA6bm[order(means_GATA6bm[,1] ),] # Order to enable a paired test
## Only Inner HMGA2, means from all cells per wells
in_wilcox_HMGA2_mw <- means_HMGA2bm %>%
  wilcox_test(HMGA2_mean ~ Setup, paired = TRUE) %>%
  adjust_pvalue(method = "BH") %>%
  add_significance("p.adj")

in_wilcox_HMGA2_mw <- in_wilcox_HMGA2_mw %>% add_xy_position()
in_wilcox_HMGA2_mw
## # A tibble: 1 × 13
##   .y.        group1      group2    n1    n2 statistic      p  p.adj p.adj.signif
##   <chr>      <chr>       <chr>  <int> <int>     <dbl>  <dbl>  <dbl> <chr>       
## 1 HMGA2_mean Above 20 µm Withi…     9     9        38 0.0742 0.0742 ns          
## # ℹ 4 more variables: y.position <dbl>, groups <named list>, xmin <dbl>,
## #   xmax <dbl>
in_co_boxes_HMGA2_mw <- ggboxplot(means_HMGA2bm, x = "Setup", y = "HMGA2_mean", 
          color = "Setup", fill = "Setup", alpha = 0.4, outlier.shape = NA) +
  xlab("Culture condition") + ylab("HMGA2 expression normalised to Hoechst signal") +
  stat_summary(fun.data = function(x) data.frame(y = 2.3, label = paste("Mean =",round(mean(x), digits = 4))), geom="text") +
 ylim(0, 3) +
  labs(title = "HMGA2 expression normalised to Hoechst signal - Inner images. Not displaying outliers", subtitle = "Paired wilcox. Labels show unadjusted p values. 1 data point =  means from all cells from 1 well") +
  scale_color_manual(values = c("#FC8D62", "#66C2A5")) +
  scale_fill_manual(values = c("#FC8D62", "#66C2A5"))

in_co_boxes_HMGA2_mw_sign <- in_co_boxes_HMGA2_mw + stat_pvalue_manual(label = "p",
                                                    in_wilcox_HMGA2_mw, tip.length = 0.01)
in_co_boxes_HMGA2_mw_sign

## GATA6

## Only Inner GATA6, means from all cells per wells
in_wilcox_GATA6_mw <- means_GATA6bm %>%
  wilcox_test(GATA6_mean ~ Setup, paired = TRUE) %>%
  adjust_pvalue(method = "BH") %>%
  add_significance("p.adj")

in_wilcox_GATA6_mw <- in_wilcox_GATA6_mw %>% add_xy_position()
in_wilcox_GATA6_mw
## # A tibble: 1 × 13
##   .y.    group1 group2    n1    n2 statistic     p p.adj p.adj.signif y.position
##   <chr>  <chr>  <chr>  <int> <int>     <dbl> <dbl> <dbl> <chr>             <dbl>
## 1 GATA6… Above… Withi…     9     9        35 0.164 0.164 ns                0.503
## # ℹ 3 more variables: groups <named list>, xmin <dbl>, xmax <dbl>
in_co_boxes_GATA6_mw <- ggboxplot(means_GATA6bm, x = "Setup", y = "GATA6_mean", 
          color = "Setup", fill = "Setup", alpha = 0.4, outlier.shape = NA) +
  xlab("Culture condition") + ylab("GATA6 expression normalised to Hoechst signal") +
  stat_summary(fun.data = function(x) data.frame(y = 1.3, label = paste("Median =",round(median(x), digits = 4))), geom="text") +
 ylim(0, 1.5) +
  labs(title = "GATA6 expression normalised to Hoechst signal - Inner images. Not displaying outliers", subtitle = "Paired wilcox. Labels show unadjusted p values. 1 data point =  means from all cells from 1 well") +
  scale_color_manual(values = c("#FC8D62", "#66C2A5")) +
  scale_fill_manual(values = c("#FC8D62", "#66C2A5"))

in_co_boxes_GATA6_mw_sign <- in_co_boxes_GATA6_mw + stat_pvalue_manual(label = "p",
                                                    in_wilcox_GATA6_mw, tip.length = 0.01)
in_co_boxes_GATA6_mw_sign

# Version2: displaying pairing datapoints from the same wells, now with the ggplot adds and not ggpubr

in_co_boxes_HMGA2_miv2 <- ggplot(means_HMGA2bm, aes(x = Setup, y = HMGA2_mean)) +
  geom_boxplot(aes(color = Setup, fill = Setup), alpha = 0.4) +
  geom_line(aes(group = Well_ID), colour = "grey") +
  geom_point(size = 1.5, aes(fill = Setup, color = Setup)) +
  xlab("Culture condition") + ylab("HMGA2 expression normalised to Hoechst signal") +
  stat_summary(fun.data = function(x) data.frame(y = 0.8, label = paste("Mean =",round(mean(x), digits = 4))), geom="text") +
  labs(title = "HMGA2 expression normalised to Hoechst signal - Inner images", subtitle = "1 data point =  means from all cells from 1 well. p values are unadjusted.") +
  scale_color_manual(values = c("#FC8D62", "#66C2A5")) +
  scale_fill_manual(values = c("#FC8D62", "#66C2A5"))

in_co_boxes_HMGA2_mi_signv2 <- in_co_boxes_HMGA2_miv2 + stat_pvalue_manual(label = "p",
                                                    in_wilcox_HMGA2_mw, tip.length = 0.01)
in_co_boxes_HMGA2_mi_signv2

## GATA6
## Only Inner GATA6, means from all cells per isolation
in_co_boxes_GATA6_miv2 <- ggplot(means_GATA6bm, aes(x = Setup, y = GATA6_mean)) +
  geom_boxplot(aes(color = Setup, fill = Setup), alpha = 0.4) +
  geom_line(aes(group = Well_ID), colour = "grey") +
  geom_point(size = 1.5, aes(fill = Setup, color = Setup)) +
  xlab("Culture condition") + ylab("GATA6 expression normalised to Hoechst signal") +
  stat_summary(fun.data = function(x) data.frame(y = 0.7, label = paste("Mean =",round(mean(x), digits = 4))), geom="text") +
  labs(title = "GATA6 expression normalised to Hoechst signal - Inner images", subtitle = "1 data point =  means from all cells from 1 well. p-values are unadjusted") +
  scale_color_manual(values = c("#FC8D62", "#66C2A5")) +
  scale_fill_manual(values = c("#FC8D62", "#66C2A5")) 

in_co_boxes_GATA6_mi_signv2 <- in_co_boxes_GATA6_miv2 + stat_pvalue_manual(label = "p",
                                                    in_wilcox_GATA6_mw, tip.length = 0.01)
in_co_boxes_GATA6_mi_signv2

Per isolation round (= mouse). OBS, only two, statistical testing is not meaningful.

## HMGA2
# Calculate means of each well

s_HMGA2og_iso <- split(HMGA2og, f = HMGA2og$Isolation)

means_HMGA2og_iso <- lapply(X = s_HMGA2og_iso, FUN = function (x) {
  split_df <- split(x, f = x$Setup)
 # head(split_df)
  lapply(X = split_df, FUN = function (a) {
    h_mean <- mean(a$HMGA2_fraction)
    h_mean
  })
}) #This function now prints out all means for each Isolation and Setup.
#means_HMGA2og
#Example for one of the wells:
means_HMGA2og_iso$iso_1
## $`Above 20 µm`
## [1] 0.7004181
## 
## $Within_20
## [1] 0.6984958
means_HMGA2bound_iso <- bind_rows(means_HMGA2og_iso, .id = "Isolation")

means_HMGA2bm_iso <- melt(means_HMGA2bound_iso, id.vars = "Isolation", value.name = "HMGA2_mean", variable.name = "Setup")

means_HMGA2bm_iso <- means_HMGA2bm_iso[order(means_HMGA2bm_iso[,1] ),] # Order to enable a paired test

## GATA6
s_GATA6og_iso <- split(GATA6og, f = GATA6og$Isolation)

means_GATA6og_iso <- lapply(X = s_GATA6og_iso, FUN = function (x) {
  split_df <- split(x, f = x$Setup)
 # head(split_df)
  lapply(X = split_df, FUN = function (a) {
    h_mean <- mean(a$GATA6_fraction)
    h_mean
  })
}) #This function now prints out all means for each Isolation and Setup.
#means_GATA6og
#Example for one of the wells:
means_GATA6og_iso$iso_1
## $`Above 20 µm`
## [1] 0.1535987
## 
## $Within_20
## [1] 0.1591578
means_GATA6bound_iso <- bind_rows(means_GATA6og_iso, .id = "Isolation")

means_GATA6bm_iso <- melt(means_GATA6bound_iso, id.vars = "Isolation", value.name = "GATA6_mean", variable.name = "Setup")

means_GATA6bm_iso <- means_GATA6bm_iso[order(means_GATA6bm_iso[,1] ),] # Order to enable a paired test
in_co_boxes_HMGA2_mi <- ggplot(means_HMGA2bm_iso, aes(x = Setup, y = HMGA2_mean)) +
  geom_boxplot(aes(color = Setup, fill = Setup), alpha = 0.4) +
  geom_line(aes(group = Isolation), colour = "grey") +
  geom_point(size = 1.5, aes(fill = Setup, color = Setup)) +
  xlab("Culture condition") + ylab("HMGA2 expression normalised to Hoechst signal") +
  stat_summary(fun.data = function(x) data.frame(y = 0.8, label = paste("Mean =",round(mean(x), digits = 4))), geom="text") +
  labs(title = "HMGA2 expression normalised to Hoechst signal - Inner images", subtitle = "1 data point =  means from all cells from 1 isolation") +
  scale_color_manual(values = c("#FC8D62", "#66C2A5")) +
  scale_fill_manual(values = c("#FC8D62", "#66C2A5"))
in_co_boxes_HMGA2_mi

## GATA6
## Only Inner GATA6, means from all cells per isolation
in_co_boxes_GATA6_mi <- ggplot(means_GATA6bm_iso, aes(x = Setup, y = GATA6_mean)) +
  geom_boxplot(aes(color = Setup, fill = Setup), alpha = 0.4) +
  geom_line(aes(group = Isolation), colour = "grey") +
  geom_point(size = 1.5, aes(fill = Setup, color = Setup)) +
  xlab("Culture condition") + ylab("GATA6 expression normalised to Hoechst signal") +
  stat_summary(fun.data = function(x) data.frame(y = 0.7, label = paste("Mean =",round(mean(x), digits = 4))), geom="text") +
  labs(title = "GATA6 expression normalised to Hoechst signal - Inner images", subtitle = "1 data point =  means from all cells from 1 isolation") +
  scale_color_manual(values = c("#FC8D62", "#66C2A5")) +
  scale_fill_manual(values = c("#FC8D62", "#66C2A5")) 

in_co_boxes_GATA6_mi

lapply(s_GATA6og, nrow)
## $iso1_well1
## [1] 4967
## 
## $iso1_well2
## [1] 7217
## 
## $iso1_well3
## [1] 11599
## 
## $iso1_well4
## [1] 8706
## 
## $iso1_well5
## [1] 8728
## 
## $iso2_well1
## [1] 1710
## 
## $iso2_well2
## [1] 2661
## 
## $iso2_well3
## [1] 1551
## 
## $iso2_well4
## [1] 1339
bind_rows(lapply(s_GATA6og, nrow))
## # A tibble: 1 × 9
##   iso1_well1 iso1_well2 iso1_well3 iso1_well4 iso1_well5 iso2_well1 iso2_well2
##        <int>      <int>      <int>      <int>      <int>      <int>      <int>
## 1       4967       7217      11599       8706       8728       1710       2661
## # ℹ 2 more variables: iso2_well3 <int>, iso2_well4 <int>
iso1 <- coculture_co[coculture_co$Isolation == "iso_1", ]

## Only Inner GATA6, grouped by the well, only from the first isolation
iso1_wilcox_gata6_g <- iso1 %>%
  group_by(Well_ID) %>%
  wilcox_test(GATA6_fraction ~ Setup) %>%
  adjust_pvalue(method = "BH") %>%
  add_significance("p.adj")

iso1_wilcox_gata6_g <- iso1_wilcox_gata6_g %>% add_xy_position()
iso1_wilcox_gata6_g[c(1:5), 11] <- 0.8
iso1_wilcox_gata6_g
## # A tibble: 5 × 14
##   Well_ID    .y.           group1 group2    n1    n2 statistic        p    p.adj
##   <chr>      <chr>         <chr>  <chr>  <int> <int>     <dbl>    <dbl>    <dbl>
## 1 iso1_well1 GATA6_fracti… Above… Withi…  4521   446    879827 8.89e- 6 2.22e- 5
## 2 iso1_well2 GATA6_fracti… Above… Withi…  6002  1215   3619262 6.84e- 1 6.84e- 1
## 3 iso1_well3 GATA6_fracti… Above… Withi…  9566  2033   9644035 5.61e- 1 6.84e- 1
## 4 iso1_well4 GATA6_fracti… Above… Withi…  7525  1181   4517786 3.55e- 1 5.92e- 1
## 5 iso1_well5 GATA6_fracti… Above… Withi…  6884  1844   5661927 1.01e-12 5.05e-12
## # ℹ 5 more variables: p.adj.signif <chr>, y.position <dbl>,
## #   groups <named list>, xmin <dbl>, xmax <dbl>
iso1_gata6 <- ggboxplot(iso1, x = "Setup", y = "GATA6_fraction", 
          color = "Setup", fill = "Setup", alpha = 0.4, outlier.shape = NA, facet.by = "Well_ID") +
  xlab("Culture condition") + ylab("GATA6 expression normalised to Hoechst signal") +
  stat_summary(fun.data = function(x) data.frame(y = 0.6, label = paste("Mean =",round(mean(x), digits = 4))), geom="text") +
 ylim(0, 1) +
  labs(title = "GATA6 expression normalised to Hoechst signal - Inner images. Not displaying outliers", subtitle = "adjusted p values per well") +
  scale_color_manual(values = c("#FC8D62", "#66C2A5")) +
  scale_fill_manual(values = c("#FC8D62", "#66C2A5"))

iso1_gata6_sign <- iso1_gata6 + stat_pvalue_manual(label = "p.adj",
                                                    iso1_wilcox_gata6_g, tip.length = 0.01)
iso1_gata6_sign
## Warning: Removed 105 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
## Warning: Removed 105 rows containing non-finite outside the scale range
## (`stat_summary()`).

## Only Inner HMGA2, grouped by the well, only from the first isolation
iso1_wilcox_hmga2_g <- iso1 %>%
  group_by(Well_ID) %>%
  wilcox_test(HMGA2_fraction ~ Setup) %>%
  adjust_pvalue(method = "BH") %>%
  add_significance("p.adj")

iso1_wilcox_hmga2_g <- iso1_wilcox_hmga2_g %>% add_xy_position()
iso1_wilcox_hmga2_g[c(1:5), 11] <- 2
iso1_wilcox_hmga2_g
## # A tibble: 5 × 14
##   Well_ID    .y.            group1 group2    n1    n2 statistic        p   p.adj
##   <chr>      <chr>          <chr>  <chr>  <int> <int>     <dbl>    <dbl>   <dbl>
## 1 iso1_well1 HMGA2_fraction Above… Withi…  4521   446   1171467 1.59e- 8 3.98e-8
## 2 iso1_well2 HMGA2_fraction Above… Withi…  6002  1215   3923418 2.85e- 5 4.75e-5
## 3 iso1_well3 HMGA2_fraction Above… Withi…  9566  2033  10583956 3.54e-10 1.77e-9
## 4 iso1_well4 HMGA2_fraction Above… Withi…  7525  1181   4419650 7.66e- 1 7.66e-1
## 5 iso1_well5 HMGA2_fraction Above… Withi…  6884  1844   6393020 6.32e- 1 7.66e-1
## # ℹ 5 more variables: p.adj.signif <chr>, y.position <dbl>,
## #   groups <named list>, xmin <dbl>, xmax <dbl>
iso1_HMGA2_g <- ggboxplot(iso1, x = "Setup", y = "HMGA2_fraction", 
          color = "Setup", fill = "Setup", alpha = 0.4, outlier.shape = NA, facet.by = "Well_ID") +
  xlab("Culture condition") + ylab("HMGA2 expression normalised to Hoechst signal") +
  stat_summary(fun.data = function(x) data.frame(y = 1.7, label = paste("Mean =",round(mean(x), digits = 4))), geom="text") +
 ylim(0, 2.3) +
  labs(title = "HMGA2 expression normalised to Hoechst signal - Inner images. Not displaying outliers", subtitle = "adjusted p values per well") +
  scale_color_manual(values = c("#FC8D62", "#66C2A5")) +
  scale_fill_manual(values = c("#FC8D62", "#66C2A5"))

iso1_HMGA2_g_sign <- iso1_HMGA2_g + stat_pvalue_manual(label = "p.adj",
                                                    iso1_wilcox_hmga2_g, tip.length = 0.01)
iso1_HMGA2_g_sign
## Warning: Removed 1470 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
## Warning: Removed 1470 rows containing non-finite outside the scale range
## (`stat_summary()`).

iso1well1 <- iso1[iso1$Well_ID == "iso1_well1", ]
iso1well2_5 <- iso1[iso1$Well_ID != "iso1_well1", ]

# For the main figure
## Only Inner GATA6,  only well1, only from the first isolation
iso1well1_wilcox_gata6_g <- iso1well1 %>%
  wilcox_test(GATA6_fraction ~ Setup) %>%
  adjust_pvalue(method = "BH") %>%
  add_significance("p.adj")

iso1well1_wilcox_gata6_g <- iso1well1_wilcox_gata6_g %>% add_xy_position()
iso1well1_wilcox_gata6_g[1, 10] <- 0.8
iso1well1_wilcox_gata6_g
## # A tibble: 1 × 13
##   .y.           group1 group2    n1    n2 statistic       p   p.adj p.adj.signif
##   <chr>         <chr>  <chr>  <int> <int>     <dbl>   <dbl>   <dbl> <chr>       
## 1 GATA6_fracti… Above… Withi…  4521   446    879827 8.89e-6 8.89e-6 ****        
## # ℹ 4 more variables: y.position <dbl>, groups <named list>, xmin <dbl>,
## #   xmax <dbl>
iso1well1$Setup <- factor(iso1well1$Setup, levels = c("Within_20", "Above 20 µm"))  #want the innermost distance group to be closest to y-axis
iso1well1_gata6 <- ggboxplot(iso1well1, x = "Setup", y = "GATA6_fraction", 
          color = "Setup", fill = "Setup", alpha = 0.4, outlier.shape = NA) +
  xlab("Distance group") + ylab("GATA6 expression normalised to Hoechst signal") +
  stat_summary(fun.data = function(x) data.frame(y = 0.6, label = paste("Mean =",round(mean(x), digits = 4))), geom="text") +
 ylim(0, 1.1) +
  labs(title = "GATA6 expression normalised to Hoechst signal - Inner images. Not displaying outliers", subtitle = "adjusted p value of well 1") +
  scale_color_manual(values = c("#66C2A5", "#FC8D62")) +
  scale_fill_manual(values = c("#66C2A5", "#FC8D62"))

iso1well1_gata6_sign <- iso1well1_gata6 + stat_pvalue_manual(label = "p.adj",
                                                    iso1well1_wilcox_gata6_g, tip.length = 0.01)
iso1well1_gata6_sign
## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_summary()`).

#As violin
violin_iso1well1_gata6 <- ggviolin(iso1well1, x = "Setup", y = "GATA6_fraction", 
          color = "Setup", fill = "Setup", alpha = 0.4, outlier.shape = NA) +
  xlab("Distance group") + ylab("GATA6 expression normalised to Hoechst signal") +
  stat_summary(fun.data = function(x) data.frame(y = 0.6, label = paste("Median =",round(median(x), digits = 4))), geom="text") +
 ylim(0, 1.9) +
  labs(title = "GATA6 expression normalised to Hoechst signal - Inner images. Not displaying outliers", subtitle = "adjusted p value of well 1") +
  scale_color_manual(values = c("#66C2A5", "#FC8D62")) +
  scale_fill_manual(values = c("#66C2A5", "#FC8D62"))

violin_iso1well1_gata6_sign <- violin_iso1well1_gata6 + stat_pvalue_manual(label = "p.adj.signif",
                                                    iso1well1_wilcox_gata6_g, tip.length = 0.01)
violin_iso1well1_gata6_sign
## Warning: Removed 53 rows containing missing values or values outside the scale range
## (`geom_violin()`).

## Only Inner HMGA2, only well1, only from the first isolation
iso1well1_wilcox_hmga2_g <- iso1well1 %>%
  wilcox_test(HMGA2_fraction ~ Setup) %>%
  adjust_pvalue(method = "BH") %>%
  add_significance("p.adj")

iso1well1_wilcox_hmga2_g <- iso1well1_wilcox_hmga2_g %>% add_xy_position()
iso1well1_wilcox_hmga2_g[1, 10] <- 1.8
iso1well1_wilcox_hmga2_g
## # A tibble: 1 × 13
##   .y.           group1 group2    n1    n2 statistic       p   p.adj p.adj.signif
##   <chr>         <chr>  <chr>  <int> <int>     <dbl>   <dbl>   <dbl> <chr>       
## 1 HMGA2_fracti… Withi… Above…   446  4521    844899 1.59e-8 1.59e-8 ****        
## # ℹ 4 more variables: y.position <dbl>, groups <named list>, xmin <dbl>,
## #   xmax <dbl>
iso1well1_HMGA2_g <- ggboxplot(iso1well1, x = "Setup", y = "HMGA2_fraction", 
          color = "Setup", fill = "Setup", alpha = 0.4, outlier.shape = NA) +
  xlab("Distance group") + ylab("HMGA2 expression normalised to Hoechst signal") +
  stat_summary(fun.data = function(x) data.frame(y = 1.6, label = paste("Mean =",round(mean(x), digits = 4))), geom="text") +
 ylim(0, 1.9) +
  labs(title = "HMGA2 expression normalised to Hoechst signal - Inner images. Not displaying outliers", subtitle = "adjusted p value of well 1") +
  scale_color_manual(values = c("#66C2A5", "#FC8D62")) +
  scale_fill_manual(values = c("#66C2A5", "#FC8D62"))

iso1well1_HMGA2_g_sign <- iso1well1_HMGA2_g + stat_pvalue_manual(label = "p.adj",
                                                    iso1well1_wilcox_hmga2_g, tip.length = 0.01)
iso1well1_HMGA2_g_sign
## Warning: Removed 50 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
## Warning: Removed 50 rows containing non-finite outside the scale range
## (`stat_summary()`).

# As violin
violin_iso1well1_HMGA2_g <- ggviolin(iso1well1, x = "Setup", y = "HMGA2_fraction", 
          color = "Setup", fill = "Setup", alpha = 0.4, outlier.shape = NA) +
  xlab("Distance group") + ylab("HMGA2 expression normalised to Hoechst signal") +
  stat_summary(fun.data = function(x) data.frame(y = 1.6, label = paste("Median =",round(median(x), digits = 4))), geom="text") +
 ylim(0, 1.9) +
  labs(title = "HMGA2 expression normalised to Hoechst signal - Inner images. Not displaying outliers", subtitle = "adjusted p value of well 1") +
  scale_color_manual(values = c("#66C2A5", "#FC8D62")) +
  scale_fill_manual(values = c("#66C2A5", "#FC8D62"))

violin_iso1well1_HMGA2_g_sign <- violin_iso1well1_HMGA2_g + stat_pvalue_manual(label = "p.adj.signif",
                                                    iso1well1_wilcox_hmga2_g, tip.length = 0.01)
violin_iso1well1_HMGA2_g_sign
## Warning: Removed 50 rows containing non-finite outside the scale range
## (`stat_ydensity()`).
## Removed 50 rows containing non-finite outside the scale range
## (`stat_summary()`).
## Warning: Removed 137 rows containing missing values or values outside the scale range
## (`geom_violin()`).

## Only Inner GATA6, wells 2-5, only from the first isolation
iso1well25_wilcox_gata6_g <- iso1well2_5 %>%
  group_by(Well_ID) %>%
  wilcox_test(GATA6_fraction ~ Setup) %>%
  adjust_pvalue(method = "BH") %>%
  add_significance("p.adj")

iso1well25_wilcox_gata6_g <- iso1well25_wilcox_gata6_g %>% add_xy_position()
iso1well25_wilcox_gata6_g[1:4, 11] <- 1
iso1well25_wilcox_gata6_g
## # A tibble: 4 × 14
##   Well_ID    .y.           group1 group2    n1    n2 statistic        p    p.adj
##   <chr>      <chr>         <chr>  <chr>  <int> <int>     <dbl>    <dbl>    <dbl>
## 1 iso1_well2 GATA6_fracti… Above… Withi…  6002  1215   3619262 6.84e- 1 6.84e- 1
## 2 iso1_well3 GATA6_fracti… Above… Withi…  9566  2033   9644035 5.61e- 1 6.84e- 1
## 3 iso1_well4 GATA6_fracti… Above… Withi…  7525  1181   4517786 3.55e- 1 6.84e- 1
## 4 iso1_well5 GATA6_fracti… Above… Withi…  6884  1844   5661927 1.01e-12 4.04e-12
## # ℹ 5 more variables: p.adj.signif <chr>, y.position <dbl>,
## #   groups <named list>, xmin <dbl>, xmax <dbl>
iso1well2_5$Setup <- factor(iso1well2_5$Setup, levels = c("Within_20", "Above 20 µm"))

iso1well25_gata6 <- ggboxplot(iso1well2_5, x = "Setup", y = "GATA6_fraction", 
          color = "Setup", fill = "Setup", alpha = 0.4, outlier.shape = NA, facet.by = "Well_ID") +
  xlab("Distance group") + ylab("GATA6 expression normalised to Hoechst signal") +
  stat_summary(fun.data = function(x) data.frame(y = 0.6, label = paste("Mean =",round(mean(x), digits = 4))), geom="text") +
 ylim(0, 2.3) +
  labs(title = "GATA6 expression normalised to Hoechst signal - Inner images. Not displaying outliers", subtitle = "adjusted p values of wells 2-5") +
  scale_color_manual(values = c("#66C2A5", "#FC8D62")) +
  scale_fill_manual(values = c("#66C2A5", "#FC8D62"))

iso1well25_gata6_sign <- iso1well25_gata6 + stat_pvalue_manual(label = "p.adj",
                                                    iso1well25_wilcox_gata6_g, tip.length = 0.01)
iso1well25_gata6_sign
## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_boxplot()`).
## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_summary()`).

#As violin
violin_iso1well25_gata6 <- ggviolin(iso1well2_5, x = "Setup", y = "GATA6_fraction", 
          color = "Setup", fill = "Setup", alpha = 0.4, outlier.shape = NA, facet.by = "Well_ID") +
  xlab("Distance group") + ylab("GATA6 expression normalised to Hoechst signal") +
  stat_summary(fun.data = function(x) data.frame(y = 0.6, label = paste("Median =",round(median(x), digits = 4))), geom="text") +
 ylim(0, 1.2) +
  labs(title = "GATA6 expression normalised to Hoechst signal - Inner images. Not displaying outliers", subtitle = "adjusted p values of wells 2-5") +
  scale_color_manual(values = c("#66C2A5", "#FC8D62")) +
  scale_fill_manual(values = c("#66C2A5", "#FC8D62"))

violin_iso1well25_gata6_sign <- violin_iso1well25_gata6 + stat_pvalue_manual(label = "p.adj.signif",
                                                    iso1well25_wilcox_gata6_g, tip.length = 0.01)
violin_iso1well25_gata6_sign
## Warning: Removed 56 rows containing non-finite outside the scale range
## (`stat_ydensity()`).
## Warning: Removed 56 rows containing non-finite outside the scale range
## (`stat_summary()`).
## Warning: Removed 215 rows containing missing values or values outside the scale range
## (`geom_violin()`).

## Only Inner HMGA2, only well1, only from the first isolation
iso1well25_wilcox_hmga2_g <- iso1well2_5 %>%
  group_by(Well_ID) %>%
  wilcox_test(HMGA2_fraction ~ Setup) %>%
  adjust_pvalue(method = "BH") %>%
  add_significance("p.adj")

iso1well25_wilcox_hmga2_g <- iso1well25_wilcox_hmga2_g %>% add_xy_position()
iso1well25_wilcox_hmga2_g[1:4, 11] <- 2
iso1well25_wilcox_hmga2_g
## # A tibble: 4 × 14
##   Well_ID    .y.            group1 group2    n1    n2 statistic        p   p.adj
##   <chr>      <chr>          <chr>  <chr>  <int> <int>     <dbl>    <dbl>   <dbl>
## 1 iso1_well2 HMGA2_fraction Withi… Above…  1215  6002   3369012 2.85e- 5 5.7 e-5
## 2 iso1_well3 HMGA2_fraction Withi… Above…  2033  9566   8863722 3.54e-10 1.42e-9
## 3 iso1_well4 HMGA2_fraction Withi… Above…  1181  7525   4467375 7.66e- 1 7.66e-1
## 4 iso1_well5 HMGA2_fraction Withi… Above…  1844  6884   6301076 6.32e- 1 7.66e-1
## # ℹ 5 more variables: p.adj.signif <chr>, y.position <dbl>,
## #   groups <named list>, xmin <dbl>, xmax <dbl>
iso1well25_HMGA2_g <- ggboxplot(iso1well2_5, x = "Setup", y = "HMGA2_fraction", 
          color = "Setup", fill = "Setup", alpha = 0.4, outlier.shape = NA, facet.by = "Well_ID") +
  xlab("Distance group") + ylab("HMGA2 expression normalised to Hoechst signal") +
  stat_summary(fun.data = function(x) data.frame(y = 1.6, label = paste("Mean =",round(mean(x), digits = 4))), geom="text") +
 ylim(0, 2.3) +
  labs(title = "HMGA2 expression normalised to Hoechst signal - Inner images. Not displaying outliers", subtitle = "adjusted p values of wells 2-5") +
  scale_color_manual(values = c("#66C2A5", "#FC8D62")) +
  scale_fill_manual(values = c("#66C2A5", "#FC8D62"))

iso1well25_HMGA2_g_sign <- iso1well25_HMGA2_g + stat_pvalue_manual(label = "p.adj",
                                                    iso1well25_wilcox_hmga2_g, tip.length = 0.01)
iso1well25_HMGA2_g_sign
## Warning: Removed 1447 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
## Warning: Removed 1447 rows containing non-finite outside the scale range
## (`stat_summary()`).

#as violin
violin_iso1well25_HMGA2_g <- ggviolin(iso1well2_5, x = "Setup", y = "HMGA2_fraction", 
          color = "Setup", fill = "Setup", alpha = 0.4, outlier.shape = NA, facet.by = "Well_ID") +
  xlab("Distance group") + ylab("HMGA2 expression normalised to Hoechst signal") +
  stat_summary(fun.data = function(x) data.frame(y = 1.6, label = paste("Median =",round(median(x), digits = 4))), geom="text") +
 ylim(0, 2.3) +
  labs(title = "HMGA2 expression normalised to Hoechst signal - Inner images. Not displaying outliers", subtitle = "adjusted p values of wells 2-5") +
  scale_color_manual(values = c("#66C2A5", "#FC8D62")) +
  scale_fill_manual(values = c("#66C2A5", "#FC8D62"))

violin_iso1well25_HMGA2_g_sign <- violin_iso1well25_HMGA2_g + stat_pvalue_manual(label = "p.adj.signif",
                                                    iso1well25_wilcox_hmga2_g, tip.length = 0.01)
violin_iso1well25_HMGA2_g_sign
## Warning: Removed 1447 rows containing non-finite outside the scale range
## (`stat_ydensity()`).
## Removed 1447 rows containing non-finite outside the scale range
## (`stat_summary()`).
## Warning: Removed 580 rows containing missing values or values outside the scale range
## (`geom_violin()`).

means_GATA6bm$isolation <- factor(nrow(means_GATA6bm))
means_GATA6bm <- means_GATA6bm %>% mutate(isolation = ifelse(grepl('iso1', means_GATA6bm$Well_ID), 'iso1', .$isolation))
means_GATA6bm <- means_GATA6bm %>% mutate(isolation = ifelse(grepl('iso2', means_GATA6bm$Well_ID), 'iso2', .$isolation))

means_gata6bm_iso1 <- means_GATA6bm[means_GATA6bm$isolation == "iso1", ]
means_gata6bm_iso1[, -4]
##       Well_ID       Setup GATA6_mean
## 1  iso1_well1 Above 20 µm 0.11492441
## 10 iso1_well1   Within_20 0.12885898
## 2  iso1_well2 Above 20 µm 0.13223626
## 11 iso1_well2   Within_20 0.13286278
## 3  iso1_well3 Above 20 µm 0.14079363
## 12 iso1_well3   Within_20 0.13628820
## 4  iso1_well4 Above 20 µm 0.10728528
## 13 iso1_well4   Within_20 0.10334277
## 5  iso1_well5 Above 20 µm 0.08952359
## 14 iso1_well5   Within_20 0.10314853
means_HMGA2bm$isolation <- factor(nrow(means_HMGA2bm))
means_HMGA2bm <- means_HMGA2bm %>% mutate(isolation = ifelse(grepl('iso1', means_HMGA2bm$Well_ID), 'iso1', .$isolation))
means_HMGA2bm <- means_HMGA2bm %>% mutate(isolation = ifelse(grepl('iso2', means_HMGA2bm$Well_ID), 'iso2', .$isolation))

means_HMGA2bm_iso1 <- means_HMGA2bm[means_HMGA2bm$isolation == "iso1", ]
means_HMGA2bm_iso1[, -4]
##       Well_ID       Setup HMGA2_mean
## 1  iso1_well1 Above 20 µm  0.4387586
## 10 iso1_well1   Within_20  0.3414553
## 2  iso1_well2 Above 20 µm  0.5468293
## 11 iso1_well2   Within_20  0.4985305
## 3  iso1_well3 Above 20 µm  0.5337545
## 12 iso1_well3   Within_20  0.4789445
## 4  iso1_well4 Above 20 µm  0.4241164
## 13 iso1_well4   Within_20  0.4239333
## 5  iso1_well5 Above 20 µm  0.4369369
## 14 iso1_well5   Within_20  0.4214000
## Only Inner GATA6, means from all cells per wells, from isolation 1
in_wilcox_GATA6_mw_iso1 <- means_gata6bm_iso1 %>%
  wilcox_test(GATA6_mean ~ Setup, paired = TRUE) %>%
  adjust_pvalue(method = "BH") %>%
  add_significance("p.adj")

in_wilcox_GATA6_mw_iso1 <- in_wilcox_GATA6_mw_iso1 %>% add_xy_position()
in_wilcox_GATA6_mw_iso1
## # A tibble: 1 × 13
##   .y.    group1 group2    n1    n2 statistic     p p.adj p.adj.signif y.position
##   <chr>  <chr>  <chr>  <int> <int>     <dbl> <dbl> <dbl> <chr>             <dbl>
## 1 GATA6… Above… Withi…     5     5         5 0.625 0.625 ns                0.142
## # ℹ 3 more variables: groups <named list>, xmin <dbl>, xmax <dbl>
## Only Inner HMGA2, means from all cells per wells from isolation 1
in_wilcox_HMGA2_mw_iso1 <- means_HMGA2bm_iso1 %>%
  wilcox_test(HMGA2_mean ~ Setup, paired = TRUE) %>%
  adjust_pvalue(method = "BH") %>%
  add_significance("p.adj")

in_wilcox_HMGA2_mw_iso1 <- in_wilcox_HMGA2_mw_iso1 %>% add_xy_position()
in_wilcox_HMGA2_mw_iso1
## # A tibble: 1 × 13
##   .y.        group1      group2    n1    n2 statistic      p  p.adj p.adj.signif
##   <chr>      <chr>       <chr>  <int> <int>     <dbl>  <dbl>  <dbl> <chr>       
## 1 HMGA2_mean Above 20 µm Withi…     5     5        15 0.0625 0.0625 ns          
## # ℹ 4 more variables: y.position <dbl>, groups <named list>, xmin <dbl>,
## #   xmax <dbl>
means_HMGA2bm_iso1$Setup <- factor(means_HMGA2bm_iso1$Setup, levels = c("Within_20", "Above 20 µm"))  #want the innermost distance group to be closest to y-axis

# Displaying pairing datapoints from the same wells, now with the ggplot adds and not ggpubr
in_co_boxes_HMGA2_miv2_iso1 <- ggplot(means_HMGA2bm_iso1, aes(x = Setup, y = HMGA2_mean)) +
  geom_boxplot(aes(color = Setup, fill = Setup), alpha = 0.4) +
  geom_line(aes(group = Well_ID), colour = "grey") +
  geom_point(size = 2, aes(fill = Setup, color = Setup)) +
  xlab("Culture condition") + ylab("HMGA2 expression normalised to Hoechst signal") +
  stat_summary(fun.data = function(x) data.frame(y = 0.65, label = paste("Mean =",round(mean(x), digits = 4))), geom="text") +
  labs(title = "HMGA2 expression normalised to Hoechst signal - Inner images", subtitle = "1 data point =  mean from all cells per well from isolation 1. p values are unadjusted.") +  
  scale_color_manual(values = c("#66C2A5", "#FC8D62")) +
  scale_fill_manual(values = c("#66C2A5", "#FC8D62")) +
  theme_clean()

in_co_boxes_HMGA2_mi_signv2_iso1 <- in_co_boxes_HMGA2_miv2_iso1 + stat_pvalue_manual(label = "p",
                                                    in_wilcox_HMGA2_mw_iso1, tip.length = 0.01)
in_co_boxes_HMGA2_mi_signv2_iso1

# As violin plot
in_co_violin_HMGA2_miv2_iso1 <- ggplot(means_HMGA2bm_iso1, aes(x = Setup, y = HMGA2_mean)) +
  geom_violin(aes(color = Setup, fill = Setup), alpha = 0.4) +
  geom_line(aes(group = Well_ID), colour = "grey") +
  geom_point(size = 2, aes(fill = Setup, color = Setup)) +
  xlab("Culture condition") + ylab("HMGA2 expression normalised to Hoechst signal") +
  stat_summary(fun.data = function(x) data.frame(y = 0.65, label = paste("Mean =",round(mean(x), digits = 4))), geom="text") +
  labs(title = "HMGA2 expression normalised to Hoechst signal - Inner images", subtitle = "1 data point =  mean from all cells per well from isolation 1. p values are unadjusted.") +  
  scale_color_manual(values = c("#66C2A5", "#FC8D62")) +
  scale_fill_manual(values = c("#66C2A5", "#FC8D62")) +
  theme_clean()

in_co_violin_HMGA2_mi_signv2_iso1 <- in_co_violin_HMGA2_miv2_iso1 + stat_pvalue_manual(label = "p",
                                                    in_wilcox_HMGA2_mw_iso1, tip.length = 0.01)
in_co_violin_HMGA2_mi_signv2_iso1

## GATA6
## Only Inner GATA6, means from all cells per isolation

means_gata6bm_iso1$Setup <- factor(means_gata6bm_iso1$Setup, levels = c("Within_20", "Above 20 µm"))  #want the innermost distance group to be closest to y-axis

in_co_boxes_GATA6_miv2_iso1 <- ggplot(means_gata6bm_iso1, aes(x = Setup, y = GATA6_mean)) +
  geom_boxplot(aes(color = Setup, fill = Setup), alpha = 0.4) +
  geom_line(aes(group = Well_ID), colour = "grey") +
  geom_point(size = 2, aes(fill = Setup, color = Setup)) +
  xlab("Culture condition") + ylab("GATA6 expression normalised to Hoechst signal") +
  stat_summary(fun.data = function(x) data.frame(y = 0.23, label = paste("Mean =",round(mean(x), digits = 4))), geom="text") +
  labs(title = "GATA6 expression normalised to Hoechst signal - Inner images", subtitle = "1 data point =  mean from all cells per well from isolation 1. p-values are unadjusted") +
  scale_color_manual(values = c("#66C2A5", "#FC8D62")) +
  scale_fill_manual(values = c("#66C2A5", "#FC8D62")) +
  theme_clean()

in_co_boxes_GATA6_mi_signv2_iso1 <- in_co_boxes_GATA6_miv2_iso1 + stat_pvalue_manual(label = "p",
                                                    in_wilcox_GATA6_mw_iso1, tip.length = 0.01)
in_co_boxes_GATA6_mi_signv2_iso1

# as violin plot
in_co_violin_GATA6_miv2_iso1 <- ggplot(means_gata6bm_iso1, aes(x = Setup, y = GATA6_mean)) +
  geom_violin(aes(color = Setup, fill = Setup), alpha = 0.4) +
  geom_line(aes(group = Well_ID), colour = "grey") +
  geom_point(size = 2, aes(fill = Setup, color = Setup)) +
  xlab("Culture condition") + ylab("GATA6 expression normalised to Hoechst signal") +
  stat_summary(fun.data = function(x) data.frame(y = 0.23, label = paste("Mean =",round(mean(x), digits = 4))), geom="text") +
  labs(title = "GATA6 expression normalised to Hoechst signal - Inner images", subtitle = "1 data point =  mean from all cells per well from isolation 1. p-values are unadjusted") +
  scale_color_manual(values = c("#66C2A5", "#FC8D62")) +
  scale_fill_manual(values = c("#66C2A5", "#FC8D62")) +
  theme_clean()

in_co_violin_GATA6_mi_signv2_iso1 <- in_co_violin_GATA6_miv2_iso1 + stat_pvalue_manual(label = "p",
                                                    in_wilcox_GATA6_mw_iso1, tip.length = 0.01)
in_co_violin_GATA6_mi_signv2_iso1

gghistogram(means_gata6bm_iso1, x = "GATA6_mean", bins = 5, color = "Setup", fill = "Setup", palette = "Set1", title = "GATA6 median per well histogram")

gghistogram(means_HMGA2bm_iso1, x = "HMGA2_mean", bins = 5, color = "Setup", fill = "Setup", palette = "Set1", title = "HMGA2 median per well histogram")

sessionInfo()
## R version 4.4.1 (2024-06-14)
## Platform: x86_64-apple-darwin20
## Running under: macOS 15.4.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: Europe/Stockholm
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] forcats_1.0.0  ggpubr_0.6.0   rstatix_0.7.2  ggplot2_3.5.2  ggthemes_5.1.0
## [6] reshape2_1.4.4 tidyr_1.3.1    dplyr_1.1.4   
## 
## loaded via a namespace (and not attached):
##  [1] sass_0.4.10        utf8_1.2.5         generics_0.1.4     stringi_1.8.7     
##  [5] digest_0.6.37      magrittr_2.0.3     evaluate_1.0.3     grid_4.4.1        
##  [9] RColorBrewer_1.1-3 fastmap_1.2.0      plyr_1.8.9         jsonlite_2.0.0    
## [13] backports_1.5.0    Formula_1.2-5      purrr_1.0.4        scales_1.4.0      
## [17] jquerylib_0.1.4    abind_1.4-8        cli_3.6.5          rlang_1.1.6       
## [21] withr_3.0.2        cachem_1.1.0       yaml_2.3.10        tools_4.4.1       
## [25] ggsignif_0.6.4     broom_1.0.8        vctrs_0.6.5        R6_2.6.1          
## [29] lifecycle_1.0.4    stringr_1.5.1      car_3.1-3          pkgconfig_2.0.3   
## [33] pillar_1.10.2      bslib_0.9.0        gtable_0.3.6       glue_1.8.0        
## [37] Rcpp_1.0.14        xfun_0.52          tibble_3.2.1       tidyselect_1.2.1  
## [41] rstudioapi_0.17.1  knitr_1.50         dichromat_2.0-0.1  farver_2.1.2      
## [45] htmltools_0.5.8.1  rmarkdown_2.29     carData_3.0-5      labeling_0.4.3    
## [49] compiler_4.4.1